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Abstract We review the orbital stability of the planar circular restricted three-body 
problem, in the case of massless particles initially located between both massive bodies. 
We present new estimates of the resonance overlap criterion and the Hill stability limit, 
and compare their predictions with detailed dynamical maps constructed with N-body 
simulations. We show that the boundary between (Hill) stable and unstable orbits is 
not smooth but characterized by a rich structure generated by the superposition of 
different mean-motion resonances which does not allow for a simple global expression 
for stability. 

We propose that, for a given perturbing mass mi and initial eccentricity e, there are 
actually two critical values of the semimajor axis. All values a < nniii are Hill-stable, 
while all values a > flungtable are unstable in the Hill sense. The first limit is given 
by the Hill-stability criterion and is a function of the eccentricity. The second limit is 
virtually insensitive to the initial eccentricity, and closely resembles a new resonance 
overlap condition (for circular orbits) developed in terms of the intersection between 
first and second-order mean-motion resonances. 
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1 Introduction 


The question of orbital stability in the circular restricted three-body problem (CR3BP) 
is a long-standing and complex problem. Our particular interest can be summarized 
in the following manner. Assume a massless particle (e.g. asteroid) orbiting a central 
star with mass toq and perturbed by a massive planet mi. We will denote by a the 
osculating semimajor axis of the particle, e its eccentricity, A the mean longitude and 
07 the longitude of perihelion. Elements with subindex 1 correspond to the perturber, 
whose orbit is considered circular (l.e. ei = 0) and exterior to that of the particle 
(a < ni). We will also assume that all motion is restricted to the plane. Under these 
considerations, given a certain eccentricity e for the particle, and fixing the angles at a 
certain value, what is the critical semimajor axis that separates the domains of stable 
and unstable motion? 

This problem has been addressed by different methods, depending on the type of 
stability under consideration. The simplest is the so-called Hill Stability, in which an 
initial condition is said to be stable if its Jacobi constant Cj is larger than the value 
Cli it acquires at the Li Euler-Lagrange point of the system. The particle will then 
be trapped within a Hill zero-velocity region that excludes the position of mi. The 
trajectory will never be able to cross the orbit of the perturber and will therefore 
remain bounded. 

The outcome of initial conditions that do not satisfy the Hill Stability criterion 
is not obvious. While the condition Cj > Cl^ is sufficient for stability, it is not 
necessary. It is possible to find solutions that do not comply with this inequality, but 
are nevertheless stable, at least for times of the order of Gyrs (e.g. Gladman 1993). As 
we will show in this paper, some of these initial conditions lie within the librational 
domain of mean-motion resonances, but others are non-resonant. 

A second estimator, this time of orbital instability, is the Resonance Overlap crite¬ 
rion, based on the work of Ghirikov (1979) and first applied to the three-body problem 
by Wisdom (1980). As its name indicates, it postulates that global chaos (and there¬ 
fore orbital instability) is triggered by the overlap of adjacent mean-motion resonances 
(MMRs). Wisdom concentrated on the case of circular orbits (e = 0) and first-order 
commensurabilities. His overlap criterion stipulates that instability will occur whenever 
the distance between two consecutive MMRs is smaller or equal to the sum of their 
libration widths. 

Using an analytical model which would later be known as the Second Eundamental 
Model for Resonance (SFMR, Henrard and Lemaitre 1983) and assuming that the sum 
of the libration widths could be approximated by twice the size of the inner separatrix 
of the commensurability farther from the perturber. Wisdom (1980) estimated that 
the averaged semimajor axis aOverlap leading to overlap may be approximated by: 


^overlap 



/ mi 

Vmo 


2/7-1 


( 1 ) 


where D = 1.3 is a constant. This expression assumes initial conditions such that the 
resonant angles a are zero, where the libration region of first-order resonances is at its 
maximum. 

Malhotra (1998) and Deck et al. (2013) presented new calculations of the overlap 
limit, using a similar analytical model description for the resonant Hamiltonian but 
with different degrees of approximations. Deck et al. (2013), for example, took into 
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account that adjacent resonances do not have the same libration half-widths, while 
both papers considered slightly different expressions for the resonance half-width. Their 
results show the same functional form in terms of the planetary mass, although with 
different numerical coefficients: D = 1.4 in the case of Malhotra (1998) and D = 1.46 
for Deck et al. (2013). 

Duncan et al. (1989) avoided analytical methods and attacked the problem using 
numerical simulations with a symplectic mapping. Their results once again yielded the 
same dependence on mi /mo, but with a larger coefficient: D = 1.49. Consequently, the 
so-called (mi/mo)^^^-law appears extremely robust to the modeling of the problem, 
although the numerical factor shows a significant spread and is less reliable. 

The case of eccentric orbits (e > 0) is more problematic. While numerical exper¬ 
iments by Quillen and Faber (2006) seem to indicate that there should not be any 
significant different in the resonant overlap limit for moderate eccentricities, analyti¬ 
cal studies of Mustill and Wyatt (2012) and Deck et al. (2013) point in the opposite 
direction. Both papers predict that even for low values of the eccentricity the overlap 
distance occurs much farther from the planets and, more surprisingly, the dependence 
with the planetary mass changes to (mi/mo)^^^. 

In this paper we revisit the resonance overlap criterion and its relation with the 
Hill stability limit. For the overlap calculations, we once again make use of the second 
fundamental model of resonance but present a new approach to the calculation of the 
overlap condition. The main difference with respect to previous studies is twofold. First, 
we show that for circular orbits there is no outer separatrix, and thus the basic idea of 
equating the resonance separation to the sum of the separatrix widths is not obvious. 

Second, we show that second-order mean-motion resonances are also important in 
determining the instability limit, even for initial conditions where their libration widths 
is minimum. Curiously, this new overlap criterion is very similar as obtained with the 
classical model, although systematically closer to the planet. We also show that our 
new overlap criterion can be used as a empirical estimate for the critical semimajor 
axis (which we denote by Ounstable) that marks the beginning of the chaotic sea and 
completely unstable orbits. 

For the Hill Stability, we present simple expressions to calculate this limit for any 
initial condition, and compare their predictions with the resonance overlap. We show 
that both criteria are different but complementary; while Cunstable marks the lower 
limit for global orbital instability. Hill stability marks the end of the stable and bounded 
motion. In between lies a rugged region of the phase space dominated by complex 
resonant structures where both stable and unstable orbits may be found. 

Finally, we analyze the case of eccentric orbits. We find that even for e ~ 0.4 
the expression for flunstable deduced for circular orbits is a very good indication of 
the global chaotic domain. The Hill Stability limit, however, is very sensitive to this 
parameter, and the transition region between both regimes grows with the eccentricity. 
Therefore a single criterion cannot be proposed as a unique law separating stable from 
unstable motion, especially for eccentric orbits. 
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2 The Resonant Hamiltonian 

2.1 Delaunay and Resonant Canonical Variables 

Since we will be working within the Hamiltonian formalism, we will introduce the usual 
modihed Delaunay canonical variables: 

L = ^/JIa ; A 

S = — \/l — e^) ; —vj (2) 

A ; Al 

where /i = Qmo and Q denotes the gravitational constant. Since the longitude of 
pericenter of the planet is constant, it does not appear as a variable of the dynamical 
system. A is the canonical momentum associated to Ai, and its value is not known a 
priori. The Hamiltonian of the system in the extended phase space can be written as: 

2 

F{L, S, A, A, 07, Al) = + mA - R{L, 5, A, w, Ai; ai, ei, 07i), (3) 

where ni is the mean motion of the perturber and R represents the disturbing function 
due to the gravitational perturbations of mi. 

Let us now suppose that the massless particle lies in the vicinity of a generic 
{p + q)/p mean-motion resonance with the perturber, such that 

(p -h q)ni - pn ~ 0, (4) 

where n is the mean motion of the particle and both p and q are positive integers. It 
is then convenient to introduce a set of resonant canonical variables (S', N, A!, a, v, Q) 
which are related to the Delaunay variables through 

S ; go- = (p -h g)Ai - pX- q-w 

N = S-L-A ■ qiy =-{p + q)\i + p\ + quji (5) 

A'=pA+{p + q)L ; qQ = X-Xi, 

where qQ is the synodic angle. The inverse transformation can be obtained easily after 
some cumbersome algebraic manipulations, and yields: 

S 

L=P(N - S) + -A' 

Q Q 

A = -^^^{N-S)-^-A' 

q q 

where M and Mi are the mean anomalies. Since the transformation (L, S, A, A, ro, Ai) - 
{S, N, A', Q) is canonical, the Hamiltonian of the extended phase space is pre¬ 
served. 


M = a + {p + q)Q 
Ml = -V +pQ 

VJ = VJx — (7 — Z/, 
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2.2 Averaging over the Synodic Angle 


In the vicinity of a mean-motion resonance, both a and v are slowly varying angles (i.e. 
long-period variables) while Q has a high frequency of the order of the orbital periods 
of the bodies. Moreover, the amplitude of the short-period variations are usually much 
smaller than their resonant and secular counterparts, and therefore have little effect 
on the long-term evolution of the system. It is thus useful to average the Hamiltonian 
with respect to Q and eliminate the short-period variations. 

The averaging is usually accomplished through a perturbation technique such as 
Hori’s method (Hori 1966, see also Ferraz-Mello 2007). Basically, we search for a Lie- 
type canonical transformation 


B{S\N*,A'\a\v\Q*) : {S,N,A',a,v,Q) {S\N*,A'*,v\Q*) 


( 7 ) 


to new (primed) variables such that the new Hamiltonian F* is independent of Q*. 
Although the construction of B is complicated when extended to high orders of the 
small parameter (here the ratio milmo), when restricted to first order it can be simply 
thought as the definite integral of F over Q in the interval [0,27r]. The same proce¬ 
dure can also be performed numerically, yielding a semi-analytical expression for the 
averaged Hamiltonian (e.g. Moons and Morbidelli 1993, Beauge 1994). 

Whichever the method adopted, we obtain a new function F*(S*, N*, A'* ^ a* ,v*) 
which is cyclic in Q*. Consequently, the corresponding canonical momenta A'* is an 
integral of motion of the system. Notice from the transformations (6) that A'* just 
appears as an additive constant in the relation between the momenta. So, independently 
of the initial conditions, we can just choose A'* = 0 and simplify both the Hamiltonian 
and the canonical transformations. 

In the averaged variables, the momenta in the Delaunay and resonant sets are 
related via: 


L* = ^{N*-S*) ; A* =- S*), 

and the averaged Hamiltonian F* can be written as: 


( 8 ) 


2p q 

( 9 ) 

where R* is now the disturbing function averaged over short-period terms. We define 
the value of L* at exact resonance as: 


/i _ (p + q) 

r * 3 “ ''^res — ni, 

res y 


( 10 ) 


and consider only initial conditions close to exact resonance. We can then expand the 
unperturbed part of F* as a Taylor series around L*res- Retaining only second-order 
terms, we can write: 

F*{S*,N*,a*,ty*) ~ -Ao{N* - S*f + Ai{N* - S*) - R*{S*, N*,a* (11) 


where Ai are positive constants that depends only on p, q and the masses: 

,2 


Ao = 


n 2 2 

2q^L* 


Ai = 


3pp 


qL* 


3 ■ 


res 


res 


( 12 ) 
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2.3 The Circular Problem 


We now consider the case where mi moves in a circular orbit (i.e. ei =0). The dis¬ 
turbing function is only function of a, the auxiliary angle v is cyclic and the associated 
momentum N* is an integral of motion. From equations (5) and (8) we can write: 

N* =S* + ^L* = = const. (13) 

P \ P J 


This implies that, given any initial values of the mean semimajor axis and eccentricity, 
their time evolution will preserve the value of N*. Both orbital elements are thus not 
independent, but coupled. The complete Hamiltonian can now be written as: 


2 2 
p q 




{p + q) 


ni{N* - S*)- R*iS*,a*-N*), 

(14) 


which is a single degree-of-freedom system parametrized by N*. 

We now turn our attention to the expression of R* (S*, a; N*) adopted for most 
analytical resonance models. From the Laplace expansion of the disturbing function, 
we will retain only the lowest-order secular and resonant terms, and thus write: 


R* = 


Qmi 

ai 


(^ffo,o(a*) + ffo,i(a*)e*^ + ffi,o(a*)e* coscr*^ . 


(15) 


In the case of first-order resonances, the expressions for the coefficients can be found 
in Brouwer and Clemence (1961) or Murray and Dermott (1999), and read: 


ffo.o(a*) 





1 

8 


2a* Doc + a*'^ dI 



(16) 




2{p + q) + a* Do 




where Doc = djda is the differential operator, and s-re Laplace coefficients. 

We will make two additional approximations. First, we will evaluate all coefficients 
at the exact resonance a^es = “res/oi- Since the perturbation is small compared to 
the unperturbed Hamiltonian, and we are only interested in the vicinity of the exact 
resonance, then the error committed here is not significant. Second, we will approximate 
the eccentricities by: 



The same arguments mentioned before are valid here, and again the error generated by 
this approximation is not relevant, at least up to eccentricities of the order of e ~ 0.5. 

Introducing these simplifications into (14), we can write the complete averaged 
resonant Hamiltonian for the circular problem as: 

F*{S*,a*;N*) = -Ao{N*-S*f+Ai{N*-S*)-CiS*-C 2 V^cosa*, (18) 
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where we have dropped constant terms, and the new coefficients are given by 


Cl = 


Qmi 

aiL*r 


90^1 (ttres) 


C 2 = 


Gmi 


Ol 






(19) 


Expression (18) constitutes a very simple analytical model for mean-motion resonances 
in the circular restricted three-body problem and, apart from the Taylor expansion of 
the unperturbed part, is identical to the Second Fundamental Model of Resonance 
(Henrard and Lemaitre 1983). 


2.4 Fixed Points and Separatrix 


For first-order (g = 1) resonances, all fixed points are located in either cr* = 0 or 
(T* = TT, and parametrized by the value of IV*. If this parameter is less than a critical 
value 


K 


1 


A1 + C1 + 



( 20 ) 


the system contains a single (stable) fixed point at a* = 0. Conversely, if N* > Nc, the 
Hamiltonian F* contains three fixed points: two centers (one located at cr* = 0 and a 
second one at a* = tt) plus one unstable point at a* = tt. The corresponding values of 
the momentum S* can be calculated analytically solving the equations of motion. The 
results, given in complex trigonometric and hyperbolic functions, can then be converted 
back to a* and e*. 

Figure 1 plots the families of fixed points for the 2/1 MMR, adopting Jupiter 
(present mass) as the perturber. The gray curves correspond to different values of 
N* = const, and the critical value N* is marked by a light black dashed curve. 
The top half-plane (positive values of e* coscr*) correspond to a* = 0, and the fixed 
points define what is usually known as the pericentric branch. All fixed points of 
the pericentric branch are linearly stable. They have low eccentricity far from exact 
resonance (shown here as a broad vertical red line), but the value of e* increases as 
a* —>■ a/es- In no case, however, does the fixed point intersect the nominal location of 
the resonance, but is always located at smaller values of the semimajor axis. 

The bottom half-plane corresponds to a* = tt. For low eccentricities, the solu¬ 
tions are again stable, and form what is known as the apocentric branch. For higher 
eccentricities (broad dashed curve) the solutions are unstable and correspond to the 
hyperbolic fixed points from which stem the separatrix of the libration regions. 

Ferraz-Mello (1988) found a simple expression relating the semimajor axis and 
eccentricity for all fixed points. Although his calculations employed the asymmetric 
expansion of the disturbing function (Ferraz-Mello 1987), the same procedure can be 
followed in the case of the SFMR. We begin writing the condition da*/dt = 0 for the 
fixed points as 


2Ao{N* - S*) - Ai - Cl 


C 2 


cos a* 


= 0 , 


( 21 ) 


where we will consider the resonant angle equal to either zero or tt. Instead of expressing 
this as an algebraic equation in S* , we recall that L* = p{N* — S*)/q from which we 
can simply obtain: 


C 2 


cos (7* 


2A0 — L* — Ai — Cl = 2A0— 

p p 



( 22 ) 
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0.6 0.61 0.62 


0.63 

* 

a /a 


0.64 0.65 0.66 


Fig. 1: Broad black curves show the families of fixed points for the 2/1 MMR in the 
(a*/ai,e*) plane, considering Jupiter as the perturber. The top half-plane (positive 
values of e* cos cr*) correspond to a* = 0, while the bottom half-plane corresponds to 
a* = TT. The position of the nominal resonant mean semimajor axis is marked by a 
broad red vertical line. The gray curves correspond to different values of N* = const 
(with A^i < . . . < N 4 ). the critical value is marked by a light black dashed curve. 


where 

J-* _ p Al -|- Cl 

" “ q 2Ao 


(23) 


constitutes the equilibrium value of the Delaunay variable. 

We next approximate \/2S* — e*\/L* ~ e* and, after some simple substi¬ 

tutions, obtain the eccentricity e* and the value of L* for all fixed points: 


-L + 

— coscr 
e* 


2Aoq 



(24) 


Values of L* < L/ give rise to the pericentric branch, while other values generate the 
apocentric and hyperbolic families. Note, however, that the value of N* is implicit 
in this equation, which may complicate the calculation of the libration width. Also, 
there is no information of the stability of each solution, which must be estimated by 
additional calculations. Finally, (24) predicts that the stable and unstable branches are 
completely symmetrical (or anti-symmetrical) with respect to L*. This is not exactly 
true, but sufficiently accurate for most purposes. 
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?L*l2L^ 

Fig. 2: Structure of the 2/1 MMR for the restricted three-body problem with Jupiter 
in a circular orbit. Thin continuous lines show the locus of stable solutions (pericentric 
branch), while broad curves mark both branches of the separatrix. The libration region 
appears only for N* > N* (shown with a dashed black curve). Note that the inner 
branch of the separatrix extends to values of e* cos cr* < 0 which implies a* = tt. 


The main advantage of (24), however, is its simplicity and ease of use. It also shows 
clearly the hyperbolic natures of the pericentric and apocentric branches and how the 
locus of fixed points tend to parabolic orbits as we approach exact resonance. 

For N* > we can calculate the borders of the libration region. These will be 
given by the values oi K = \/2S* cos a*, with a* = 0, tt such that the Hamiltonian 
coincides with the value at the hyperbolic hxed point. Together with the value of N* 
we can then transform them into orbital elements and calculate the values: (a*n,e*jj) 
and (oJutieout)- The hrst pair will define the branch of the separatrix separating the 
libration zone from the inner circulation region or, more correctly, its intersection with 
the a* = 0, TT axis, while the second will mark the appearance of the outer circulation 
domain. 

Figure 2 shows the structure of the MMR with Jupiter as the perturber (present 
mass). The red curves were calculated using the SFMR Hamiltonian (18), while the 
black curves show results using a semi-analytical model for the resonant Hamiltonian 
in which the averaged disturbing function is evaluated numerically at every point. 
The SFMR shows very good agreement with the exact calculation, especially for low 
eccentricities, which justifies the use of the analytical model for near-circular orbits. 
However, the SFMR systematically underestimates the libration width, a datum that 
will be taken into consideration in later stages of this work. 
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To simplify our notation, we will denote by inner separatrix the locus of points 
separating the inner circulation region from the librational domain, calculated for all 
N* > Nc- Similarly, we will refer to the boundary between the libration and outer 
circulation domains as the outer separatrix of the resonance. Note that the inner 
branch of the separatrix extends below e* cos cr* = 0, indicating that it is also present 
for a* = tt; thus, the limit of the librational domain for circular orbits (e* = 0) is not 
given by the curve N* = N* but for larger values of the integral of motion. 

An important feature of Figure 2 that has received little attention is the fact that 
the outer separatrix does not extend down to circular orbits. Although for every value 
of N* > Nc the system contains both an inner and outer branch, these appear in the 
(o*,e*) plane with different eccentricities. This is a natural outcome of the fact that 
both sets (a*n, e*^) and (Ooutj fiout) must lead to the same value of N*. 

If we analyze this figure for a given (fixed) value of the eccentricity, then for e* > 
0.18 both branches of the separatrix will be present and the librational domain will be 
contained within. For e* < 0.18, however, there is no outer branch of the separatrix. In 
other words, for this range of eccentricities all values of the semimajor axis a* < Ores 
will yield N* < Nc and will thus correspond to a circulation. There is in fact a libration 
region, but is much smaller, and limited by the inner separatrix and the value of the 
semimajor axis such that N* = Nc- 


3 The Classical Resonance Overlap Criterion 

Wisdom (1980) established what may be called the condition for the classical resonance 
overlap in the CR3BP. Following the ideas of Walker and Ford (1969) and Chirikov 
(1979), overlap is said to occur when, for a given value of the mean eccentricity e*, the 
distance Aa*cs between adjacent mean-motions resonances is smaller or equal to the 
sum of their libration widths ZiUgep. The same idea was also adopted by later analytical 
calculations, such as Malhotra (1998) and Deck et al. (2013), although this last paper 
improved the model taking into account that the half-widths between neighboring 
MMRs are not equal. 

All these estimates also employed the SFMR as the resonance model, although with 
different degrees of approximation of the Hamiltonian, or in the deduction of the reso¬ 
nance width. In many cases these approximations were made in order to obtain results 
in simple explicit expressions, even though the errors introduced were not adequately 
checked. 

Finally, resonance overlap was calculated in the representative plane defined by 
a* = zu = 0 and evaluated at e* = 0. These choices of initial conditions also allowed 
to consider only the role of first-order commensurabilities. Second-order MMRs have a 
minimum libration size for a* = 0, while the separatrix width of third-order resonances 
tends to zero for circular orbits. So, focusing on first-order resonances appears a good 
approximation, especially when studying overlap for circular orbits. 

In this section we once again deduce the resonance overlap limit for circular orbits 
(i.e. e* =0), following the same assumptions as described above. We will, however, 
reduce the approximations to a minimum, even if this implies results which are not in 
explicit analytical expressions. 
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3.1 Libration Width for First-Order Resonances 

For any given first-order mean-motion commensurability (i.e. q = 1), we wish the 
estimate the width of the libration domain for circular orbits. In other words, we are 
interested in the distance Z\aJep(p) between the nominal position of the (p -|- l)/p 
resonance and the edge of the inner branch of the separatrix at e* = 0. 

We begin with the expression (18) for the Hamiltonian of the SFMR. Since the 
value of N* is constant, we can add a quantity equal to C\N* with no change in the 
dynamics. Thus, we obtain a new expression given by: 

F* = -Ao{N* - S*f + {Ai + Ci){N* - S*) - coscr* (25) 

or, writing {N* - S*) = L*/p, 

F* = + -C2V^cosa*. (26) 

p^ p 

We have kept the denomination of this function, although Hamiltonians (25) and (26) 
differ by a constant amount. We can simplify this expression even further. Using the 
dehnition of LJ given by (23) we can write {Ai -|- Ci) = ‘lA^L^jp. Introducing this 
equality in the Hamiltonian, and after dividing both sides by a factor Aq jp , we obtain 

F* = ^F* =-L*^+ 2L* L*-^^^V^cosa*. (27) 

Aq Ao 

To calculate the libration width for a given value of N* we must first determine the 

value of the Hamiltonian at the hyperbolic hxed point. From (24) we can write the 

value of S’^ and L’^ at a given fixed point as: 

\ (28) 

where both S’^ and are function of N*. Substituting into (27), the Hamiltonian 
of the hyperbolic fixed point will be given by: 

i^hyper(A^*) = -i1v' + 2L:L^ + ^y^. (29) 

We now search for those values of L* = Lggp and S* = such that the 

Hamiltonian at cr* = 0 attains the same value. This is simply: 

= LgV" + 2Lh l:,p - ^(30) 

Note, however, that TJep a-nd 6'sep are not independent, but related through the chosen 
value of N*. Without any loss of generality, we can then write: 

L:,p = L% + AL* ; s:,p = Sh + AS = Sh - , (31) 

where the last equality stems from the constraint that N* = const. Since we are 
interested in calculating the libration half-width for circular orbits (i.e. S'sep = 0), the 
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second equation reduces to S'jy = AL*/p. Introducing these relations into expressions 
(29) and (30) and equating the values of both Hamiltonians, we obtain 

{AL*f + 2{L*n - L*) AL* + V2Ai? = 0, (32) 

which admits the non-trivial solution: 


AL* = l:,^-L*^=P^ 


36*2 

4Aoy 


(33) 


We can now calculate the distance AL^ep between the inner separatrix and the 
exact resonance as 


ALl,^ = LLp - (LLp - L*m) + {L*m - Ll) = (^1 + (34) 

From (12) and (19) we can write 


36*2 


J-'c ^ 6 * \ 

r> 2 ^res 5l,o(<^resj 
mo 2p^ 


0.4 mi 

p mo 


j- *3/2 

-^C 5 


(35) 


where we have used the approximation 5i,o(ci:?es) — ~0.8p (Malhotra 1998). Last 
of all, converting the results to semimajor axis, we obtain 


Aatepip) — 1-6 ai 


mo) 


2/3 / 


mo 


,mo + mi 


1/3 


(p+ 1)2/3’ 


(36) 


which is very similar to the expression given by Wisdom (1980). However, in order to 
preserve accuracy even for low-degree resonances, we will avoid any analytical trans¬ 
formation between p its value of a*es- 


3.2 Separation Between First-Order Resonances 


For this step we will follow the deduction presented in Morbidelli (1999) although, once 
again, trying to circumvent any non-essential simplifications that may affect accuracy 
of our model. 

The nominal location of the {p -|- l)/p and (p -|- 2)/(p -|- 1) MMRs can be written 
as; 


^res 


mo 


mo + mi 


1/3 


p + 1 


2/3 


from which their separation is given by: 


a 


*1 

res 


\mo + miJ \p + 2J 

(37) 


Aa*es{p) = ai^ 


mo 


mo + mi 


1/3 


p + 1 

LV^ 


2/3 



(38) 
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We now assume that p ^ 1, and expand each of the terms inside the square brackets 
as a Taylor series. Retaining only first-order terms, we can approximately write: 


P 

p + 1 

p + 1 
p + 2 


2/3 

2/3 


^ x2/3^^ 2 1 

P+1/ ~ 3(p+l) 

1 1 
p + 2/ 3 (p + 2) 


(39) 


Finally, introducing both expressions into (38), and after some simple algebra, we 
obtain: 


(p) 


2^ / mo 1 

3^ V™o + mi/ (p-h l)(p + 2) ■ 


(40) 


This is moderately different from the expression in Wisdom (1980) and Deck et al. 
(2013) which adopt 2\a*es = (2cti)/(3(p-|-1)^)- On one hand, it is important to explic¬ 
itly keep the mass ratio between the perturber and the star, which might be important 
for massive planets, or even binary stellar systems. On the other hand, the difference in 
the dependence on p is also significant, especially for low-degree resonances. It is simple 
to see that while the approximate relation adopted by Wisdom (1980) systematically 
overestimates the true separation, the opposite occurs for (40). We found that a much 
more accurate estimate, even for low values of p may be written as: 


^ares(p) 


2a f mo 1 

3^^\mo + mi) (p-h l)(p-h 3/2) ’ 


(41) 


in other words, half way between both of the original approximations. 


3.3 The Classical Overlap Condition 

Having expressions for the libration half-width (eq. (36)) and the separation between 
consecutive first-order resonances (eq. (41)), we can now proceed to calculate the con¬ 
dition for resonance overlap. Following Wisdom (1980), this is said to occur whenever 
the order of the resonance acquires a value p = Pc such that: 

^ares(Pc) = 2Z\a*ep(pc). (42) 

Deck et al. (2013) improved this estimation replacing the right-hand side with Z\a/ep(Pc) + 
Z\o/ep(Pc + 1), leading to a significant change in the results. Nevertheless, as will be 
discussed further on, this modification of the classical condition will not prove necessary 
and we will adopt expression (42). 

Whatever our choice, the idea remains the same: overlap is defined when the inner 
separatrix of the (pc -|- l)/pc MMR intersects the outer branch of the (pc -|- 2)/(pc -|- 
1) resonance. As in the original Chirikov criteria, the crossing of separatrices of two 
commensurabilities generates a dynamical route leading to orbital instability. However, 
it is important the stress that mere proximity between resonances is not sufficient; 
overlap requires the existence and intersection of two separatrix. If both are not present, 
then overlap cannot be said to occur. 
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Solving equation (42), we obtain the value of pc that signals resonance overlap (for 
circular orbits) according to this model as: 


4.8pc(Pc + 3/2)(pe + l)^/^ 



(43) 


Since giving results in values of pc is awkward, we can transform them to values of a* 
using expression (37). As before, we avoid analytical simplihcations and any a-priori 
assumption that pc ^ 1, but choose to perform a one-dimensional numerical fit of the 
critical semimajor axis as function of mil mo. The result yields: 


^overlap — 


1 - 1.225 ( — 

mo 


0 . 28 - 


(44) 


This, then, is the minimum initial mean semimajor axis, for circular orbits, such that 
the massless particle lies in an unstable region generated by the overlap of first-order 
interior MMRs with a perturbing planet (also in circular orbit) of mass mi. 

It is important to stress that this calculation, as well as the resulting expression, 
is given in averaged orbital elements. This is vital when attempting to compare its 
predictions with numerical simulations of the exact Newtonian differential equations. 
To obtain the overlap limit in osculating semimajor axis we must perform the back 
transformation, which yields the following approximate relation 


^overlap — ai 


1 - 1.06 



(45) 


This expression can now be compared directly with numerical integrations, as long as 
the initial angles are chosen equal to zero. Note that the exponents of both descriptions 
of the overlap criterion are similar to 2/7 ~ 0.286, but slightly smaller. 


4 Numerical Tests of the Overlap Criteria 

4.1 The Ae Indicator 

An interesting tool to analyze the behavior of planetary systems is the so-called Max¬ 
imum Eccentricity Method (MEM) (e.g. Dvorak et al. 2004), which follows the maxi¬ 
mum value of e attained by a given body during a numerical simulation. For initially 
eccentric orbits, a better indicator is the difference between the maximum and mini¬ 
mum values, or the value of Ae. 

Although it is not a measure of chaotic motion, this Ae indicator is an extremely 
useful tool to map the resonant structure in N-body problems. Its application is very 
simple. A grid of initial conditions in a representative plane is integrated numerically 
for a timespan T larger than the period of the slowest angle of the system. During the 
integration we keep track of the minimum and maximum values attained by one of the 
actions (call it J), and calculate the difference A.J = Jmax ~ •/min- Finally, we plot 
the value of AJ as a color graph in the plane of initial conditions. 

Figure 3 shows an example for the simple pendulum. The left-hand frame shows 
the level curves of the Hamiltonian in the (J, 9) plane. Since the system has only a 
single degree of freedom, the plane is the complete phase space. The stable fixed point 
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Fig. 3: Left: Phase plane of the simple pendulum. Right: Maximum increment of the 
momentum J for a grid of initial conditions in the plane. 


appears at the center of the graph {J = 0, 6 = 0), while the separatrix is highlighted 
in dashed lines. 

We now divide the plane in a 100 X 100 grid of initial values of the action and the 
angle. Each initial condition is then integrated for a total time T. The resulting values 
of J are shown in a color plot on the right-hand side of the figure. Initial conditions 
close to the stable fixed point appear with a dark color, indicating very small values of 
AJ. Conversely, initial conditions close to the separatrix appear as very light colored 
regions, corresponding to large changes in the momentum. Thus, even if we were not 
able to obtain explicitly the phase curves of the Hamiltonian (i.e. plot on the left), the 
color plot on the right would allow us to estimate both the location of possible fixed 
points as well as the separatrix of the resonance region. 

Having defined the method, we can apply it to a mean-motion resonance in the 
circular restricted three-body problem. Figure 4 shows the resulting map for the 2/1 
MMR in the representative plane {a/ai,e) of (osculating) initial conditions, where all 
angles were chosen initially equal to zero. All initial conditions were integrated with a 
Bulirsch-Stoer based N-body code for T = 10® orbits of the perturber, and the color 
code is the difference in eccentricity attained by each particle (i.e. Zie). Thin gray lines 
mark different values of N*, with Nc highlighted with a dashed curve. 

We can now compare these results to the predictions of the semi-analytical reso¬ 
nance model. The predicted pericentric branch is shown as a dashed black curve, while 
both separatrix branches are indicated with broad black lines. The agreement is very 
good, indicating that the semi-analytical model is a very precise tool for estimating the 
features of a given mean-motion resonance. More importantly, it shows that, for circu¬ 
lar orbits, there is no outer separatrix for the 2/1 MMRs, and the same holds for any 
other first-order commensurability. For e = 0, all initial conditions with a < flres show 
practically no change in the eccentricity and remain close to circular orbits throughout 
the integration. Only those initial conditions with a > Ores show an increase in the 
eccentricity, reaching a maximum at the location of the inner branch of the separatrix, 
roughly located at a ~ 0.64ai. 
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Fig. 4: Values of Ae in the representative plane (a/ai,e) of osculating initial condi¬ 
tions in the vicinity of the 2/1 MMR with Jupiter (current mass). All angular variables 
were initially chosen equal to zero. The pericentric branch of zero-amplitude solutions 
is shown as a dashed black curve and both branches of the separatrix in broad black 
continuous lines. Curves of constant N* are indicated with gray lines, with Nc high¬ 
lighted as a dashed curve. The regions of maximum variation of the eccentricity appear 
red, while those associated with small changes are indicated in blue. 


Again, the nonexistence of the outer separatrix for low eccentricities is the result of 
plotting the representative plane in non-canonical variables. For all values of A^* > Nc 
both branches of the separatrix exist, although with different values of the semimajor 
axis and eccentricity. 


4.2 Large-scale Dynamical Map of the Representative Plane 

In order to compare the different versions of the resonance overlap criterion, we con¬ 
structed a dynamical map in the plane (o/oi, mi/mo), with 800 values of the osculat¬ 
ing semimajor axis and 200 values of the mass ratio between mi/mo G [10~®, 10~^]. 
Initial osculating eccentricities and angles where chosen equal to zero, and all orbits 
were integrated for T = 10^ orbital periods of the perturber. Results are shown in Fig¬ 
ure 5, where blue regions mark small changes in the eccentricity, and dark red indicate 
hyperbolic motion and/or escapes. 
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Fig. 5: Ae map of the {a/ai,mi/mo) plane of initial conditions (in osculating ele¬ 
ments) with circular orbits. Final hyperbolic orbits are shown in dark red. The colored 
curves are different predictions of resonance overlap: Wisdom’s criterion with D = 1.30 
(black), the analytical expression by Deck et al. (2013) with D = 1.46 (blue), numeri¬ 
cal fit by Duncan et al. (1989) with D = 1.49 (red), and the criterion developed here 
(cyan). The location of the main first-order MMRs are indicated with white text. 


The map shows the different resonances present in this interval of semimajor axis, 
starting with the 2/1 MMR at a/ai — 0.63, up to first-order resonances with p —I 20 
for semimajor axis ratios above 0.95. These resonances generate a “saw-tooth” shape 
limit for the unstable region. As we will see further on, not all correspond to first-order 
commensurabilities. The same figure shows, in continuous curves, the predictions of the 
different overlap criteria discussed in the previous section. To allow for an adequate 
comparison, all the values of a*.^,gj.iap were transformed to osculating semimajor axes. 

While all the criteria appear similar, the expressions deduced by Duncan et al. 
(1989) and Deck et al. (2013) seem a good fit to the lower extrema of the saw-tooth 
structure. In the case of the value given by Duncan et al. (1989), this is not unexpected 
since it was obtained by a least-square fit of a series of numerical simulations using a 
symplectic map. The black curve, corresponding to the original prediction by Wisdom 
(1980), passes through the middle of the saw-tooth region, as though these structures 
were evened out. The value D = 1.30 then appears to yield a better “average” boundary 
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between the stable and unstable domains. Last of all, the broad cyan curve presents 
the predictions of our model, as given by eq. (45). Although it has a different functional 
form, it yields practically the same results as Wisdom (1980) for, say, mi/wo < 10”^. 
For larger masses, however, there is an increasing and noticeable discrepancy, and our 
values are systematically closer to the perturber. 

Before pursuing a more detailed comparison between these models, it is important 
to review whether the assumptions behind the resonance overlap criterion are in fact 
consistent with the structures shown in the dynamical maps. First, there is evidence 
of second-order resonances (5/3, 7/5, 9/7, etc.), particularly close to the instability 
limit, that interact with the first-order MMRs and contribute to generate the chaotic 
domain. These are not considered in the overlap criterion developed so far. 

A second feature that can be perceived from the dynamical map is that the outer 
separatrix of the first-order MMRs appears to have little effect in the dynamics of the 
system. The reason behind this was mentioned before: the outer separatrix (located 
at a < Ores) only exists for eccentricities above a certain threshold, and is not present 
for circular orbits. Then, the only truly resonant region for circular orbits occurs for 
a > Ures and is characterized by the inner separatrix. Thus, the idea of using the 
libration width for the inner separatrix as a proxy for the outer branch is questionable, 
at least for circular orbits. 


4.3 New Criterion with C* and 2"'^-order MMRs 

To analyze this point in more detail, Figure 6 shows the same dynamical map as before, 
only this time we have superimposed the nominal position of the first and second- 
order MMRs (dashed black and green lines, respectively) in terms of the osculating 
semimajor axis. While a/es is a weak function of the perturbing mass (and actually 
decreases for larger values of mi), the difference between a* and a shows a much 
stronger dependence. Moreover, our choice of angles for the representative plane (i.e. 
A = Al = 07 = 0) corresponds to the maximum value attained by the semimajor axis 
due to short-period oscillations, which leads to exact resonance occurring closer to the 
planet for increasing mi. 

The same figure also shows, in continuous black curves, the inner separatrix for the 
first-order commensurabilities. We can now see a much clearer agreement between the 
resonant structure and the results of the numerical simulations. As we consider larger 
values of mi, the region close to the inner separatrix of each first-order resonance 
becomes increasingly chaotic, generating a region of orbital instability, characterized 
by Ae —1. Concurrently, a different chaotic domain appears, linked to second-order 
commensurabilities, whose effects have been almost negligible until that point. The 
size of both chaotic regions increase with the perturbing mass until, at some point, 
both intersect. Global chaos then sets in and most initial conditions between both 
commensurabilities become unstable. 

A complementary view is presented in Figure 7, in the form of new maps in the 
(a/ai,e) plane for five values of mi/mo. As before, dark regions are associated to 
small changes in the eccentricity during the total integration time, while the opposite 
occurs for initial conditions identified with lighter tones. The location of second-order 
resonances are indicated in gray lines, while the separatrix of first-order commensura¬ 
bilities are shown in color curves. These were calculated using a semi-analytical model 
(e.g. Beauge 1994) and not with the SFMR to allow for a better correlation with 
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Fig. 6: Same dynamical map shown in the previous figures, this time superposed with 
the resonant structure of first (black) and second-order (green) MMRs. The location 
of each pericentric branch in shown in dashed lines, while the extension of the inner 
separatrix for circular orbits (first-order resonances only) is shown in continuous lines. 
The Cyan curve shows the result of the new resonance overlap criterion developed using 
these resonances. The broad white curve, corresponding to eq. (49), is the approximate 
location of the beginning of the global unstable region. 


the numerical results. It is important to keep in mind, however, that the separatrix 
have been drawn assuming isolated resonances and, as such, do not take into account 
perturbations form nearby commensurabilities. 

While the dynamical features are complex and show evidence of the interaction 
between many different resonances, some characteristics may be deciphered analyzing 
their evolution as function of the planetary mass. For mi/mo = 10”^'^, all initially 
circular orbits with ajax > 0.88 are unstable, a value close to the intersection point 
between the 11/9 resonance and the inner separatrix of the 5/4 commensurability. 
Initial conditions with smaller semimajor axis appear primarily regular, although some 
chaotic motion is visible In the vicinity of a/ai ~ 0.85, associated to a near-intersection 
between the 9/7 resonance and the inner separatrix of the 4/3. 

Between mi/mo = I0~^'^ and mi/mo = 10~^'^,the 9/7 and 4/3 MMRs intersect 
and all circular initial conditions closer to the planet are unstable. A new isolated region 
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Fig. 7: Top left-hand frame shows a zoom of the dynamical map in the (o/ai, mijmo) 
plane presented in the previous figure. Five values of the perturbing mass are high¬ 
lighted with horizontal black lines. The intersection between adjacent first and second- 
order resonances are indicated with black circles. The remaining graphs present new 
dynamical maps in the (a/ai,e) plane for each value of mi/mo. The location of the 
second-order resonances are indicated with gray lines, while the separatrix of first-order 
MMRs are shown in different colors. Each commensurability is indicated on the top of 
the frames. 


of chaotic motion appears at a/ai ~ 0.80, generated by the interaction between the 7/5 
commensurability and the inner separatrix of the 3/2 resonance. For mi/mo = 10~^ ® 
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this chaotic layer becomes more extensive and finally merges with the global instability 
region for mi/mo — 10”^'^. 

The same behavior can be observed for other values of the planetary mass and 
semimajor axis. Thus, it appear that global instability may be estimated by the inter¬ 
section of the inner separatrix (i.e. located at o > Ores) of a given first-order MMR with 
the nominal location of its adjacent second-order commensurability (i.e. black circles 
in the top left-hand frame of Figure 7). In fact, the libration width of the second-order 
resonance does not appear important at all. This result enormously simplifies the cal¬ 
culations, since we do not require to model the resonance structure of the second-order 
resonances; all we require is to estimate its location in the semimajor axis domain for 
any given planetary mass. 

These numerical results led us to postulate a new overlap criterion based in the 
interaction between adjacent first and second-order MMRs. To calculate this new crite¬ 
rion, we first write average mean motion of a generic first-order MMR and its adjacent 
second-order commensurability as: 


p + 1 


ni 


2p 3 
2p+l 


ni. 


(46) 


This expression excludes “false” second-order resonances such as the 4/2 or 6/4, since 
both the numerator and denominator present in n*res are odd numbers. The distance 
between them can be calculated using the same procedure as in the Section 3.2, which 
now gives: 


^Ores(Pc) = ai ^ 


mo 


mo + mi 


1/3 


2pc + i y^^ 

2pc + 3 / 


Pc 


Pc + 1 


2/3- 


(47) 


For the libration half-width of the first-order commensurabilities we use expression 
(36), analogous to our derivation of the classical overlap criterion. Equating (47) with 
(36) we can solve for pc such that Aa*{pc) = Z\a/es(Pc)- After transforming the result 
in terms of the osculating semimajor axis a, we finally obtain our new overlap condition 


^overlap — 0-1 


1-0.91 


mi 

mo 


(48) 


where both numerical coefficient were determined from a non-linear regression with val¬ 
ues of mi/mo G [10~^, 10”^]. Figure 6 shows, with a broad Cyan curve, the prediction 
of this new criterion. Since the calculations are not exact, there is a slight difference 
with respect to the actual resonance intersections. It is interesting to note that this 
result is very similar to the one determined using the classical overlap condition (i.e. 
Figure 5), although both were obtained with completely different assumptions. 

Since the SFMR underestimates the separatrix width (see Figure 2), the predictions 
of expression (48) fall short of the instability limit as shown in Figure 7. However, we 
found that it is possible to improve the estimation simply changing the numerical 
factor. The white curve in Figure 6 corresponds to the equation 


^unstable — 


1 - 0.75 


mi 

mo 


0.26n 


(49) 


which has the same functional dependence with the perturbing mass as the new overlap 
criterion, but with a smaller coefficient. Although this is an empirical correction of our 
analytical result, it results in a fairly accurate estimate for the value for the osculating 
semimajor axis that marks the beginning of the region of global instability. 
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5 The Hill Stability Criterion 

A different stability criterion may be defined in terms of the value of the Jacobi constant 
Cj of the particle, as compared with its value Cl^ at the Li Lagrange point. If 
Cj > Cli, then the massless body is forever trapped in a zero-velocity curve that 
encompasses the central mass toq, but which does not contain the perturbing mass 
mi. In such a case the initial condition is said to be Hill Stable. Conversely, if the 
initial conditions are such that Cj < Cl^, then the stability is not guaranteed, and 
may suffer close approaches with mi and become temporarily trapped by the smaller 
primary. Its mo-centric motion will then be characterized by hyperbolic orbits and 
thereby cataloged as unstable. 

The Jacobi constant then constitutes a valuable asset with which to determine the 
(Hill) stability of a given massless particle in the realm of the CR3BP. We will make 
use of this feature to develop a complementary stability criterion that will later be 
compared with the predictions of the resonance overlap criterion. 


5.1 Hill Stability in Orbital Elements 


The main problem with the practical application of this criterion is to calculate both the 
Jacobi constant for Li and its value for any given initial condition in orbital elements. 
Equations usually found in the literature either make use of series expansions (which 
may or not yield accurate results for large planetary masses) or require Cartesian 
coordinates and velocities in the rotating reference frame. In the next paragraphs we 
will introduce a way to overcome both limitations. 

Seidov (2004) introduced closed formulas to calculate the position and Jacobi con¬ 
stant for Li. Let us call the distance from mo to Li. We can then write 

I'Ll = ai(l - 5 li) (50) 


where the new auxiliary quantity will satisfy 5li —> 0 as mi —?■ 0. Seidov (2004) found 
that Sli and the planetary mass are related through: 

mi _ (1 - (5li)^(1 + Sli + 

mo“ 5li3-3dL,+Sl^) ■ ^ > 

Although this equation is exact for any value of mi, it is implicit in Therefore, 
in order to determine the position of Li as function of the planetary mass, we must 
solve it using successive approximations. Since Sl^ is usually a small quantity, we can 
rewrite (51) is a form more adequate for the iterative process. This is: 










(52) 


Choosing zero as the initial guess in the right-hand side, this expression can be solved 
in just a few iterations, and yields very precise results even for large values of the 
perturbing mass. 

Having determined the location of the Lagrange point, the value of the Jacobi 
constant Cl^ can also be calculated using a closed formula. Again following Seidov 
(2004), we can write: 


Cl 


g{mo + mi) / 3 - 12e -b 15e^ - lOe^ - 4e^ ^ 


m 


(l-2e-e2)2 


(53) 
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where e = 5 li(1 - Sl^)- 

Our next task is the obtain an explicit equation for Cj for any massless particle 
in terms of its orbital elements. Although the Jacobi constant is usually expressed 
in Cartesian coordinates in the rotating reference frame, it is also possible to write 
it in an inertial frame (e.g. Murray and Dermott 1999). Assuming initial conditions 
A = Ai = ro = 0 (i.e. the particle is located in its pericenter and all bodies are aligned), 
then we find: 



where ni is the mean motion of the perturbing planet, and Ai are the distances from 
the massless body to rrii, given by: 

Zii = ai — a(l — e). 

With these simple expressions we may now calculate the boundary of the Hill 
stability limit. Given values of mo, mi and e, we search for the value of the semimajor 
axis flHill such that Cj = Cl^- This procedure can be performed for a range of 
planetary masses mi, from which we can construct our stability limit curve in the 
(a/ai,mi/mo) plane. 

An important point is that, contrary to the resonance overlap criterion, all the 
expressions developed here are given in terms of osculating orbital elements. The only 
drawback is that this criterion cannot be expressed in a simple closed form for arbitrary 
values of mi. Although some approximations may be found in the literature (e.g. 
Gladman 1993, Deck et al. 2013), usually based on the Hill description of the restricted 
three-body problem, they are usually only accurate for small perturbing masses (of 
the order of mi/mo ~ 10”^ and lower), and even then just in the circular case. In 
comparison, the model developed above is still purely analytical and extremely fast to 
use once implemented in a computer code. 


5.2 Application to Circular Orbits 

In order to compare the predictions of this criterion with those obtained from resonance 
overlap, we will first analyze the case e = 0. Results are shown in Figure 8, again in 
the foreground of the same dynamical map discussed in the previous section. However, 
we have also highlighted in gray all initial conditions that are unstable within the 
integration timespan. These include cases where the final orbit is hyperbolic as well as 
those that cross the orbit of the perturber while retaining e < I. The broad white curve 
is the empirical proxy for Uunstable, while the broad black curve is the prediction of 
the present Hill Stability Criterion flHill • All versions of the resonance overlap criterion 
are shown in thin continuous lines. 

Figure 9 shows a zoom of the region around the spike near the 3/2 MMR. This plot 
was constructed from a fresh set of numerical simulations covering a 300x300 grid in 
semimajor axis and mass ratios. As before, all angles were initially set equal to zero 
and the integration time for each initial condition was T = 10^ orbits of the perturber. 
Each spike shows a complex shape on the left side, but a significantly more smooth 
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Fig. 8: Dynamical map for initially circular orbits (l.e. e = 0) and all angular variables 
equal to zero, for a grid of initial conditions in the (a/ai, mi/mo) plane. This map 
is equivalent to the one discussed extensively in the previous section, except that all 
unstable initial conditions are highlighted in gray. The thin curves mark the predictions 
of the different resonance overlap criteria; Black: Wisdom (1980), Blue: Deck et al. 
(2013) and Red: Duncan et al. (1989), Cyan: new overlap limit defined in section 4.3. 
The broad black curve shows the predictions of the Hill Stability Criterion, while the 
broad white curve is the empirical global instability limit Ounstable(™i/™o) given by 
eq. (49). 


edge on the right side. The origin of the dichotomy is not clear, although it could be 
caused by higher-order resonances that are only visible in this level of detail. 

As always, thin continuous lines shows the result of the different resonance overlap 
criteria, while the broad white curve is the proxy for the beginning of the global Insta¬ 
bility region. Except for a big wedge around mi/mo ~ 10~ ' , this prediction of the 
beginning of the chaotic sea appears very accurate. Again, with the exception of the 
wedge, all initial conditions with a > flunstable (he. beyond the broad white curve) are 
unstable. 

The broad black curve shows the results of the Hill Stability limit. Although the 
saw-tooth structure seems to be indifferent to the presence of the limit curve, the values 
of Ae are different on both sides. All orbits satisfying the condition a > aniii are in 
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Fig. 9: Same as previous Figure, but now zooming in on a spike centered roughly around 
the 3/2 mean-motion resonance. Notice the change in stability on both sides of the Hill 
limit. 


fact unstable, lending credibility to our expressions in terms of the orbital elements. 
For smaller semimajor axes, however, the increase in eccentricity is significantly lower, 
indicating that, although there is a certain excitation of the orbit, it is bounded and 
not sufficient to render the motion unstable. 

It is interesting to note that the region located between both stability criteria 
(i.e. ttHiii < a < aunstable) IS also characterized by sharp borders between stable 
and unstable orbits. More importantly, the condition Cj > Cl^ is sufficient but not 
necessary for stability; initial conditions exist for which motion is always bounded to 
the central mass even though the zero-velocity surface includes both primaries. 


5.3 Application to Eccentric Orbits 

Figure 10 shows the results of the calculation of the Hill Stability limit, in the (a/oi, mi / mo) 
plane, for different values of e. The boundary for low masses is very sensitive to the 
initial eccentricity, while the change seems less drastic for larger masses. Even so, even 
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Fig. 10: The Hill Stability limit for different values of the particle’s eccentricity (iden¬ 
tified by the number accompanying each curve). 


a low eccentricity introduces a significant change in the stability limit, which should 
be noticeable in a dynamical map. 

Figure 11 now shows four dynamical maps of Ae where, contrary to Figure 8, 
we have not highlighted orbits that are unstable while maintaining low eccentricities. 
Each frame corresponds to a different initial eccentricity of the particle, as indicated 
in white text. As the eccentricity increases, so does the instability region close to 
the planet, even for small perturbing masses. The libration width of each resonance 
also increases, generating islands of very regular motion immersed in regions of highly 
chaotic behavior. Many high-degree MMRs are noticeable for a —> ai, as well as for 
semimajor axis below that corresponding to the 2/1 commensurability. The resonance 
structure thus becomes increasingly complex, much more so than observed for circular 
orbits. 

The broad black curve indicates the Hill Stability limit, as calculated for each 
value of the eccentricity. We can see that it follows very closely the inner edge of the 
chaotic domain, and it appears fairly clear that all initial conditions with a < auiii 
are stable, at least for the total integration time (again set at T = 10® orbits of mi). 
The broad Cyan curves correspond to the empirical value of flunstable developed in this 
work, as calculated for circular orbits. Surprisingly, the same functional form shows a 
very good agreement with the beginning of the global chaotic region, independently of 
the eccentricity of the particle. Thus, it appears that this upper stability limit is an 
extremely weak function of e, even up to e = 0.4. 
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Fig. 11: Dynamical maps of Ae in the (a/ai,mi/mo) plane of initial conditions for 
four different values of the particles’ eccentricities (indicated within each frame). Small 
changes in the eccentricity after 10® years integration time are shown in blue, while 
unstable orbits leading to ejection are shown in red. In each plot, flunstable in shown 
with a broad Cyan curve, while the Hill Stability limit (calculated for the specific value 
of e) is shown in black. The dashed blue line shows the analytical derivation of the 
overlap criterion for eccentric orbits developed by Deck et al. (2013). 


Finally, the dashed blue line shows the predictions of the overlap criterion developed 
for eccentric orbits by Mustill and Wyatt (2012) and later refined by Deck et al. (2013). 
Using a simple analytical model for the growth of the libration domain as function 
of the eccentricity, both papers deduce that the overlap limit has a functional form 
proportional to (e(7Tii/mo))^^®, although they differ in the value of the numerical 
constant. Here we have adopted the version of Deck et al. (2013). 
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Fig. 12: Dynamical map of Z\e for a grid of initial conditions in the (a/ai, e) plane of 
osculating elements (all angles equal to zero), with Jupiter as the perturber. Integration 
time was set to T = 10^ orbits of the perturber. The broad continuous curves show 
the Hill stability limit (black), the overlap limit deduced in this paper (cyan) and the 
value of Ounstable (white). The dashed green curve is the prediction of the eccentric 
overlap criterion by Deck et al. (2013). 


Trying to extend the overlap criterion to non-circular orbits is complex. As we 
showed previously, the outer separatrix of first-order resonance does not exist for e ~ 0 
but does appear for some minimum value, and thus should be taken into consideration 
at some point. Third-order resonances also become noticeable for eccentric orbits and 
begin to overlap with lower-order commensurabilities for sufficiently high eccentricities. 

Since it is difficult to take into consideration all these factors, it is not surprising 
that the analytical model does not show a good agreement with the dynamical maps. 
In fact, they do appear to be an “average” between both the circular resonance overlap 
and Hill stability limits. Curiously, however, for eccentricities e < 0.1 there is a very 
good coincidence with the Hill Stability curve. 

Figure 12 shows a different example. This time we constructed a Zie dynamical map 
in the (a/ai,e) representative plane, considering Jupiter (current mass and circular 
orbit) as the perturber. We defined a grid of initial conditions with 2048 values of 
semimajor axis and 576 values of the initial eccentricity. All angles were taken equal to 
zero. The total integration time was only T = 10^ orbits of the perturber, not sufficient 
to eject many unstable orbits, but sufficient to show large increases in the eccentricity. 
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Fig. 13: Numerical integration of a particle with a = 0.8, e = 0.4 and all angles initially 
set to zero. The perturbing mass was chosen equal to mi/mo = 10~®. The particle 
exhibits a random walk between adjacent first and second-order MMRs before finally 
exceeding the perturber’s orbit at f ~ 2.5 X 10^ years. During the complete integration 
the eccentricity remains close to its initial value with Ae of the order of 0.1. 


Apart from the main first-order commensurabilities (2/1, 3/2 and 4/3), there is dis¬ 
tinct evidence of second and third-order resonances, particularly the 7/4 at a/ai ~ 0.69 
and the 8/5 at a/ai ~ 0.73. These are negligible for near-circular orbits, but become 
increasingly important for higher eccentricities (e.g. Bodman and Quillen 2014). At 
e ~ 0.2 these begin to overlap and generate a large chaotic region close to the Hill 
Stability curve, shown here again as a broad black curve. 
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Fig. 14: Similar to Figure 11, but now the abscissa is a function of the semimajor 
ajcis ratio defined by: u{a/ai) = —[(a/ai)^/^ — 1)]”^. In this new variable all the 
first-order resonances appear equidistant, thus allowing a better visualization of the 
resonant structure closer to the perturber. As before, the black curve shows the Hill 
stability limit, while cyan corresponds to Ounstable- The eccentric overlap criterion of 
Deck et al. (2013) is shown as a dashed blue curve. 


Gladman (1993) also reported, from the result of a few numerical simulations, that 
some orbits above the Hill stability curve may survive for long time-spans as long as 
they satisfy the Hill condition for circular orbits. As we show here, the relationship 
between instability and the Hill limit is more complex. 

The broad cyan curve shows our version of the resonance overlap criterion, as 
deduced for circular orbits, while the white vertical line shows the value of Uunstable- 
In particular, the value of Ooverlap shows a very good agreement with the beginning of 
the global chaotic sea, and appears fairly independent of the eccentricity. Finally, the 
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dashed green curve corresponds to the eccentric overlap criterion as proposed by Deck 
et al. (2013). 

Figures 11 and 12 show that the region between nniii and Ounstable grows with 
increasing value of e and lower values of mi. This transition region is characterized 
by a complex resonant structure and shows the existence of both stable and unstable 
orbits. Stable motion is usually found deep inside the libration domain of mean-motion 
resonances, while unstable orbits abound elsewhere. 

For perturbing masses below mi/mo ~ 10^ — 10^, the results of Figure 11 seem 
to indicate that particles with a > aniii are not very excited and remain with low- 
to-moderate eccentricities. To investigate this in more detail, Figure 13 shows the 
dynamical evolution of a single initial condition for a total integration time of 3 X 10^ 
years. The perturbing mass was chosen equal to mi/mo = 10~®. As mentioned in 
the caption, the body suffers several jumps in semimajor axis, becoming temporarily 
trapped between adjacent first and second-order mean-motion resonances. After ap¬ 
proximately 25000 years the orbit finally surpasses the semimajor axis of the perturber 
and eventually diffuses to the outer regions of the system. The eccentricity, however, 
remains bounded throughout the integration, with a maximum increase (with respect 
to its initial value) of the order of Ae — 0.1. 

Since orbital instability is not necessarily associated to hyperbolic orbits, Figure 
14 repeats the results of Figure 11, where the gray region now shows all orbits that 
became unstable within T = 10^ orbits, either reaching e > 1 or becoming planet 
crossers. We also plotted the data in terms of an auxiliary function u{a/ai) instead of 
the semimajor axis. This function is defined by: 


u{a/ai) 



(56) 


In this new variable all first-order resonances appear equidistant, thus allowing us to 
have a better visualization of the structure of the representative plane closer to the 
perturbing mass. 

As before, the lower bound of the chaotic sea is very well represented by the circular 
estimate of Ounstable, while the Hill Stability correctly delimits the appearance of the 
resonance forest. The eccentric overlap formula of Deck et al. (2013) is seen in dashed 
blue lines. As before, this later criterion shows a good agreement with the Hill limit for 
low eccentricities and planetary masses above mi/mo ~ 10”^, but rapidly diverges 
for both lower masses and higher eccentricities. 


6 Conclusions 

In this paper we presented the results of a series of high-resolution Ae dynamical maps 
in a representative plane of the planar circular restricted three-body problem, whose 
aims where two-fold: (i) obtain a detailed visualization of the limit between stable and 
unstable orbits (in the Hill sense) and (ii) estimate the resonant structure of mean- 
motion commensurabilities and their relation with the stability limit. These results 
were used to test the predictions of the Hill stability limit and different versions of the 
resonance overlap criterion (circular and elliptic case). 

In particular, using the Second Fundamental Model for Resonance (Henrard and 
Lemaitre 1983), we obtained an alternative derivation of the overlap criterion based 
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on the interaction between first and second-order commensurabilities which appears 
to be responsible for the transition between local and global chaotic motion. The re¬ 
sulting expression is reminiscent of the one presented by Wisdom (1980), although 
with a smaller numerical coefficient and slightly different functional dependence on the 
perturbing mass. 

Chirikov’s postulation of the overlap criterion was constructed for systems in which 
the interacting resonances occurred in the same set of canonical variables and, implic¬ 
itly, with the same expressions for the integrals of motion. These conditions do not 
apply for the problem at hand where each mean-motion resonance is characterized by 
distinct expressions of the integral N* and, consequently, different set of canonical 
variables. Classical studies have circumvented this problem by avoiding canonical vari¬ 
ables and analyzing their interaction in the (a, e) plane, variables which are common 
to all MMRs. However, as we was seen throughout this paper, these orbital elements 
are not adequate and the outer separatrix does not appear for quasi-circular motion, 
leading to a series of assumptions that have not always been justihed. Thus, although 
Chirikov’s overlap criterion has been widely used in the circular restricted three-body 
problem, its applicability has not always been adequately verified. 

It is important to stress that the aim of this paper has not been to introduce a 
new rigorous theory for resonance overlap, but to present a qualitative criteria that 
preserves the principles of the classical formulation but (more importantly) reproduces 
the dynamical features observed in the numerical simulations. 

Finally, we showed that it is not possible to characterize the stability limit using a 
single criterion. The Hill Stability criterion defines the maximum value of the semimajor 
axis (for a given eccentricity) for which all initial conditions are (Hill) stable. On the 
other hand, an empirical variation of our resonance overlap limit (49) is a useful proxy 
to estimate the minimum semimajor axis for which all initial conditions are unstable. 
In between lies a transition domain characterized by many resonant islands and rich 
in both stable and unstable motion. The size of this region grows with increasing 
eccentricity or lower values of the perturbing mass, and may occupy a signihcant portion 
of the phase space. 
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